Chapter 10 

An Analytical Application: Optimization 
of a Stirling Engine Based on the Schmidt 
Analysis and on the Adiabatic Analysis 


Alles soli so einfach wie moglich 
gemacht werden, aber nicht einfacher. 

Everything should be made 
as simple as possible, but not simpler. 

Albert Einstein 


10.1 Introduction 

Stirling engines are external combustion engines converting thermal energy into 
mechanical energy by alternately compressing and expanding a fixed quantity of air 
or other gas (called the working or operating fluid) at different temperatures [126]. 
Stirling engines were invented by Robert and James Stirling in 1818. Despite their 
high efficiency and quiet operation they have not imposed themselves over the Diesel 
and Otto engines. In recent years interest in Stirling engines has grown, since they are 
good candidates to become the core component of micro Combined Heat and Power 
(CHP) units. In this chapter, we discuss an optimization experiment performed on 
Stirling engines. In particular, optimization algorithms are applied to the Schmidt 
and to the adiabatic analyses. These are two simple and rather idealized analytical 
models of the Stirling machine. Before discussing the optimization issue we briefly 
recall the basic elements of the Stirling cycle, and the Schmidt and the adiabatic 
analyses. 
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Fig. 10.1 The ideal Stirling cycle 


10.1.1 The Stirling Thermodynamic Cycle 


Stirling engines are based on the Stirling regenerative thermodynamic cycle which 
is composed of four thermodynamic transformations: 

• an isothermal expansion at high temperature, 

• an isochoric regenerative heat removal, 

• an isothermal compression at low temperature, 

• an isochoric regenerative heat addition. 

Since the operating fluid is expanded at high temperature and compressed at low 
temperature a net conversion of heat into work is attained. The theoretical efficiency 
of the cycle in case of complete reversibility equals that of the ideal Carnot cycle, 
as stated by the Reitlinger theorem [127]. An ideal Stirling cycle between the tem¬ 
peratures 7/ and Th (7) < 7^), and between the volumes V/ and Vh (V/ < 14 ) is 
represented in Fig. 10.1 and is described by the following equations 


W U2 = Ji pdV = MRT h In ^ > 0 

W 2>3 = o 

W 3>4 = MRT, In \ < 0 


21.2 = Wl ,2 = MRTh In ijy >« 

22.3 = Mc v (7/ - T h ) < 0 

23.4 = W 3 ,4 = MRT, In < 0 


W 4 J =0 24,1 = Mc v (T h - 7/) = -22,3 > 0 

Wne, = W U2 - W 3> 4 = MR (T h - 7}) In ^ > 0 r, = ^ = 1 - ^ = Vcarnot 

( 10 . 1 ) 

where W m , n and Q m ,n respectively are the amount of work and the heat exchanged 
by the system during the transformation from the status m to the status n, p is the 
pressure, and M the mass of the operating fluid in the system, R is the specific 
gas constant, c v is the specific heat at constant volume of the gas, W net is the net 
work output, and rj the thermodynamic efficiency of the cycle. 02,3 and Q\,\ are 
exchanged regeneratively, thus they are not included into the efficiency equation. 
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Cooler Regenerator Healer 


Fig. 10.2 Stirling engine schematic representation 


10.1.2 The Schmidt Analysis 

Schmidt analysis [128-130] is an ideal isothermal nonlinear model for the simulation 
of Stirling machines. 

The working space of a Stirling machine is composed of: 

• a compression space (c), 

• a cooler (k), 

• a regenerator (r), 

• a heater ( h ), 

• an expansion space ( e ). 

Figure 10.2 is a schematic representation of a Stirling machine and its spaces and 
pistons. The fluid flows back and forth between the expansion and the compression 
spaces crossing the heater first, then the regenerator, and finally the cooler. The fluid 
is displaced by the motion of a piston (the displacer) and is compressed and expanded 
by the motion of another piston (the power piston). 

The main assumptions of the Schmidt analysis are: 

• constant thermodynamic properties of the operating fluid, 

• sinusoidal volume variations in the expansion and the compression spaces due to 
the pistons motion, 

K (0) = V d ,e + ^ (1 + COS 0), Vc (0) = V d , c + ^ (1 + cos (0 - a)) 

(10.2) 

• constant volume of the heater, the regenerator, and the cooler, 

• constant and uniform temperature equal to T d in the expansion space and in the 
heater, 

• constant and uniform temperature equal to Tk in the compression space and in the 
cooler, 

• constant and linearly varying temperature in the regenerator between Tk and Th , 

• uniform pressure in the whole working space, 
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p ( 6) = MR 



+ 



Th ~ Tk 



+ 



(10.3) 


0 E [0, 27r] defines the actual phase in the cycle, a the phase lag between the volume 
variation in the expansion and in the compression space. V stands for the volume, 
Vd for the dead volume, V sw for the swept volume, p for the pressure, M for the total 
mass of operating fluid, R for the specific gas constant, T for the thermodynamic 
temperature. The subscripts e and c stand for expansion and compression spaces 
respectively, the subscripts h,r,k for the heater, the regenerator, and the cooler. The 
regenerator mean effective temperature is defined 


T r = 


Th ~ Tk 

In ^ 
111 T k 


(10.4) 


The following nondimensional parameters are used in the analysis: 

• the temperature ratio r = ^ 

• the volume ratio ip = 


v. 


sw,c 


Vs 


sw,e 


the regenerator dead volume ratio x r = 
the hot dead volume ratio Xh = Vd f ^ Vh , 


V r 


sw,e 


sw,e 


the cold dead volume ratio Xk = Vd y c + Vk . 


sw,e 


Substituting Eqs. 10.2 into 10.3 yields 


MRTk ip r 

——-= (1 + COS (6> - a)) + - (1 + cos 0) + H 

P (0) V sw ,e 2 2 


(10.5) 


where 


H = 


x r r In - 

_ _ T_ 

1 — T 


+ XhT + Xk 


( 10 . 6 ) 


is the reduced dead volume. The phase 
in the cycle is such that 

tan Oq — 


angle 6>o at which the pressure is minimum 

ip since 

—T7-• (10.7) 

r + y; cos a 


Defining 


K = 


2MRT k 

Vsw,e 


Y = T + ip + 2H 


( 10 . 8 ) 


Equation 10.5 can be written in the form 


p(d) = 


K 

Z(1 +5 cos (6»-6> 0 )) 


(10.9) 
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where 

7r2 + cos a 

r + ip + 2H 

is the pressure swing ratio. The mean pressure over the cycle is 


1 /' 27r K 

Pm = - / p (0) d0 = — . 

F 2tt 7o 

The expansion and compression work during one cycle is given by 
w r> Js A]/ P™ dVg TTPm T JU , e^sin^o 

W e = Qe=<ppdV e = p —d6= - 

/ Jo ^ 1 + Vl -s 2 

... „ f ... I' 27 " dV c np m ipv sw , e 5 sin (<9 0 - a) 

W c =Q c =<pp dV c = p —d6= - - 

J J o d9 1 + Vl — 5 2 


It follows that the net power output and the efficiency of the cycle are 


( 10 . 10 ) 


( 10 . 11 ) 


( 10 . 12 ) 

(10.13) 


W net = w e + w c = 


nPmVsw,ed(l ~ T) sin 6>p 

i + Vi -s 2 


T)=l-T = T] carn ot • (10.14) 


Thus, the Schmidt analysis still yields the ideal Carnot efficiency. The work output 
depends upon the following parameters: x r , Xh, Xk, 'ip, t , a, M , R, Tk, V sw , e . 

W net can be expressed in nondimensional form by dividing by MRT\ or by 

Pmax Vsw,e 


W net = 


w, 


net 


2tt5 (1 — r) sin Oq 


MRTk V1 — S 2 (t + 'ip + 2^1 + Vl — 


Wnet = 


w, 


nS \/1 — S (1 — r) sin Oq 
Pmax Vs w, tot V1 + 5^1 T f ip\ ^1 T V1 — 


(10.15) 


(10.16) 


where V sw ,tot = Vsw,e + Vsw,c- The net work output given by the nondimensional 
Schmidt analysis just depends upon x r , Xh, Xk, V t, a. According to the Schmidt 
analysis, the dead volumes are always reducing the work output, and the smaller is 
r the higher is the net work output. Thus, for a given r value, the optimal values of 
the parameters are x r = Xh = Xk = 0. From this situation it follows that a meaning¬ 
ful nondimensional optimization would involve just two input variables: ip and a. 
However, it must be considered that the optimum configurations also depend upon 
x r , Xh, x k , and r, since W net = W net (r, ip, a , H), and W net = W net (r, ip, a, H). 
In fact, all the terms in Eqs. 10.15 and 10.16 can be written as functions of r, ip, a. 
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77, and 77 depends upon x r ,Xh,Xk. Figure 10.3 shows the nondimensional net work 
output as a function of ip and a. 

In real engines, actually, there is no point in removing the regenerator and the 
heat exchangers, since even if their volume is a “dead” volume, not being swept by 
the pistons, their presence is fundamental for the engine to work properly. In fact, 
being an external combustion engine, the heat exchangers are the only thermal energy 
source and sink. 


10.1.3 The Adiabatic Analysis 

Schmidt’s hypothesis that the expansion and the compression spaces are assumed to 
be isothermal, as a consequence of the cycle being reversible, implies that the heat 
is exchanged directly by these spaces with the two sources. The regenerator is also 
ideal. Therefore, all the heat transfer processes occurring in the real world, do not 
influence the Schmidt analysis. 

The adiabatic analysis is a sort of improved Schmidt analysis where the expansion 
and compression spaces are assumed to be adiabatic. In this way, the heat enters 
and leaves the engine only through the heat exchangers which are distinct from 
the expansion and compression spaces. The adiabatic analysis is still an idealized 
nonlinear model of Stirling engines since it retains the assumption of ideal (i.e. 
reversible) heat exchangers and regenerator. This is still quite a heavy assumption 
since the heat exchangers and the regenerator are the core of Stirling machines. 
Therefore the adiabatic analysis can still give quite erroneous results, even if more 
realistic than Schmidt analysis, and predicts an overall engine efficiency not too far 
from the one of the Carnot cycle. 

The adiabatic assumption makes it impossible to obtain a closed form solution, 
as it was for the Schmidt analysis and demands an iterative solving procedure to be 
enforced. 

The main assumptions of the adiabatic analysis are: 

• the thermodynamic properties of the operating fluid are constant, 

• the engine consists of five spaces: the expansion space ( e ), the heater ( h ), the 
regenerator (r), the cooler (7), the compression space (c) (See Fig. 10.4), 

• the volume variations in the expansion and compression spaces are sinusoidal and 
follow Eq. 10.2, 

• the volumes in the heater ( 14 ), the regenerator (V r ), and the cooler ( 14 ) are 
constant, 

• the temperatures in the heater (7^), and the cooler (7^) are constant and uniform, 

• the temperature in the regenerator is constant and linearly varying between 7^ 
and Th , thus, the regenerator mean effective temperature is given by Eq. 10.4, 

• the expansion and the compression spaces are adiabatic, 

• the pressure is uniform within the working space (p = p e = ph = p r = Pk = Pc ), 
and, under the ideal gas equation, is expressed as 
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(a) Wnet , contour view 



(c) Wn e t , contour view 



(b) Wnet , 3d view 



(d) Wnet , 3d view 


Fig. 10.3 Nondimensional net work output according to the Schmidt analysis as a function of tJj 
and a, for r = \,x h = = ^ (# = 0.188) 
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Fig. 10.4 The Stirling engine working space 


p (9) = MR 


Vc(0) 
T c (0 ) 


+ 





vyyy' 

T e (0)) 


(10.17) 


where 6 is the phase angle, M the overall mass of the operating fluid, R the 
thermodynamic constant, specific of the operating gas. 

Solving the adiabatic analysis means to compute, for each value of the crank angle 
6 , the volume, the temperature, and the mass of operating fluid for each engine section, 
and the pressure in the working space. The amount of heat and work exchanged during 
the cycle are finally computed. 

From the above assumptions, the adiabatic analysis depends upon ten variables: 
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• the volume and the temperature of the expansion and the compression spaces 
(V e , V c , T e , T c \ 

• the mass of the operating fluid in each space (M c , M k , M r , Mh, M e ), 

• the pressure in the working space (p). 

Thus, in order to solve the adiabatic analysis, ten equations are needed, and they are: 

• two volume variation equations for the expansion and the compression spaces 
(Eq. 10.2), 

• two energy balance equations for the expansion and the compression spaces, 

• five state equations, one for each space, 

• one continuity equation. 

For solving the energy balance equations we need to compute the mass flow rates in 
and out of the expansion and the compression spaces. We designate M' e ^ h the mass 
flow rate from the expansion space to the heater, and M k ^ c the mass flow rate from 
the cooler to the compression space. We also define the upwind temperature at the 
interfaces T e ^ k , T k ^ c which are conditional on the direction of the flow 


T e if M' e ^ h > 0 

Th ifM^<0 


T k if M' k ^ c > 0 
T c if M[^c < 0. 


(10.18) 


The state equation for a generic space and the continuity equation in differential form 
can be written as 


dp 

P 


+ 


dV 

V 


d M t d T 

~W + Y" 


(10.19) 


d M e + d M h + d M r + d M k + d M c = 0 (10.20) 


respectively, where d M e = —M f e ^ h and d M c = M' k ^ . The energy equations for a 
generic space is 

d Q + c p T in M' in - c p T out M' out = dW + c v d(MT) (10.21) 

which, for the expansion and the compression spaces becomes 

-c p T e ^ h M' e ^ h = pdV e + c v d(M e Te) (10.22) 

c p T k ^ c M' k ^ c = pdV c + c v d(M c T c ). (10.23) 

Here, d Q and dW stand for infinitely small quantities of transferred heat and work, 
c v and c p are the specific heats of the operating fluid at constant volume and at 
constant pressure respectively. With a few algebraic passages Eqs. 10.22 and 10.23 
can be written in the forms 
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d M c 


p dV c + p 

RTk^c 


(10.25) 


Q 

where 7 = R = c p — c v . By differentiating the state equation the following are 
derived for the heater, the regenerator, and the cooler, respectively: 


M h V h Mr V r M k V k 

d M h = - dp = —dp, d M r = —dp = —dp, dM k = —dp = —dp 

P RT h p RT r p RTk 

(10.26) 

Substituting Eqs. 10.24-10.26 into the continuity equation, with a few algebraic 
passages, yields 



(10.27) 


From the state equation, the following equations hold for the expansion and the 
compression spaces: 



dVe 

V e 




d Vc 
V c 



(10.28) 


Applying the energy equation to the heater, the cooler, and the regenerator it is 
possible to express the amount of heat exchanged by each section: 

d Q h = C -P^-?- - Cp (T e ^ h M' e ^ h - T h M' h ^ r ) (10.29) 

d Qr = -- Cp {T h M' h ^ r - T k M' r ^ k ) (10.30) 

d Q k = C -^^~ - c p {T k M' r ^ k - TWM'^). (10.31) 

K 

Finally, the work done is given by 

dW e = pdV e , dW c = pdV c , dW ne t=dW e +dW c . (10.32) 


The choice of the operating fluid determines R , c p , c v , 7. A crank angle step size 
A 9 must be defined. The steps of the adiabatic analysis, from iteration n to iteration 
n + 1 , are: 


update the crank angle 0 ^+!) — Q( n ) _j_ A 0 , 

update the values of the expansion and the compression volumes (Eq. 10.2) and 


their derivatives V e ( ' !+1) , V c ( " +1) , (^f) ( \ (^) 


(n+D 
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• update 7+ +1) , 7+ +1) , wj" +1) , W c ( " +1) , Q 
integration, 

• update the conditional temperatures T^ h \ 


O+i) nO+l) /)0+l) 
h ’ 


7’A+ I ' (Eq. 10.18), 


by numerical 


• update the pressure p (/7+1) (Eq. 10.17) and the pressure derivative 



O+i) 


(Eq. 10.27), 

• for each engine space, update the mass M^ nJrV> (using the ideal gas equation), the 

mass derivative )^ +1 ^ (Eqs. 10.24-10.26), and the mass flow M' according to 
the following equation 


M' e ^ h = -dM e , M' h ^ r = -dM e - dM h , M' r ^ k = dM c + dM k , M' k ^ c =dM ( 


update 


(We\ 

V d ^ ) 

(d &V" +1) 


(n+D 


/dW c \ 

V / 


O+l) 


\ 0 + 1 ) 

(d T c\ 

) •' 

I" j 


0 + 1 ) 


(10.33) 


(Eq. 10.28), 


/ d Wnet \ 

9 V do ) 


O+l) 


( d Qh \ 

V / 


0+D 


(dQr\ 

V d ^ ) 


O+l) 


, and 


\ d6 ) 


(Eqs. 10.29-10.32) 


The adiabatic analysis, even though it is not an initial value problem, is solved as 
an initial value problem and the procedure is started by inputting a set of arbitrary 
conditions at t = 0 (or 0 = 0) and fixing a At (or AO) step. The process above 
is repeated up to convergence, that is, up to when the initial transitory has been 
dampened out. Experience has shown that the most sensitive measure of convergence 
is the residual regenerator heat Q r at the end of the cycle, which should be zero. 
The first order explicit Euler method [131] yields fairly accurate results for Stirling 
adiabatic analysis. For better accuracy, it is suggested to employ the fourth order 
Runge-Kutta method [132]. 


10.2 The Case 

Typically, the optimization of any thermal engine aims at 

• the maximization of the power output P ou t of the engine, 

• the maximization of the thermodynamic efficiency r] of the engine. 

Here, we address Stirling engines design by means of Schmidt and adiabatic analyses 
with these two objectives in mind. Although Schmidt and adiabatic simulations 
are idealized models, they are a good and cheap starting point for Stirling engines 
evaluation. 

A Stirling engine is said to be of /3 type when the power piston is arranged 
within the same cylinder and on the same shaft as the displacer piston. Two other 
engine configurations are possible: a when two power pistons are arranged in two 
separated cylinders, and 7 when a power piston and a displacer are arranged in 
separate cylinders. The type of the Stirling engine does not influence the Schmidt 
and the adiabatic analyses. Figure 10.5 shows a very simple model of a Stirling 
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Fig. 10.5 j3 Stirling engine with rhombic drive 


engine of /3 type with a rhombic drive. The rhombic drive is one of the most popular 
drive mechanisms for Stirling engines. In the figure are shown (top-down): 

• the heater (in dark red), 

• the hot cylinder wall (in pink), 

• the regenerator (in yellow), 

• the displacer piston (in green), 

• the cooler (in blue), 

• the cold cylinder wall (in grey), 

• the power piston (in orange), 

• the rhombic drive (in light purple), 

• the drive gears (in cyan). 

The idea behind the optimization experiment performed is to find the optimal 
configurations according to the two analytical models, and to compare the differences 
in the results. Of course, before running the simulations we need to define a few 
constraints on the engine. These constraints are necessary otherwise the optimization 
process would result in engines with, for instance, infinite hot temperature, volumes, 
and power output. The basic assumptions for the exercise are: 

• Schmidt and adiabatic analysis are accepted as valid means for Stirling engine 
simulations and are employed in the optimization process. This is equivalent to 
accept the assumptions discussed in Sects. 10.1.2, 10.1.3, 

• the objective of the optimization are the maximization of the engine’s power output 
and thermodynamic efficiency 



rpm 

— W net s r\ ’ 

60 


W net 

Qin 


(10.34) 


where rpm stands for the engine frequency in revolutions per minute, Q\ n is the 
heat input to the engine in one cycle, W ne t the work output given by the engine in 
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one cycle. The simulation code employed to perform the Schmidt and the adiabatic 
analyses was written using C++ language. 


10.3 Methodological Aspects 

In this section we discuss briefly and in chronological order the choices made for the 
setup of the optimization process. We roughly retrace the same steps seen in Chap. 8. 
Most of the observations made in that chapter are still valid for this application, and 
will not be repeated here. A schematic representation of the choices which have been 
made, and which will be discussed, is given in Fig. 10.6. 


10.3.1 Experiments Versus Simulations 

As usual there are two possible ways for collecting data: by means of experiments, 
or by means of simulations. In the case of a generic and introductory approach for 
sizing Stirling engines, there is no way to adopt the experimental approach: in a real 
engine, in fact, a very large number of parameters come into play, and we need some 
other means for collecting data quickly. 

As anticipated in Sect. 10.1, we choose to address the optimization of Stirling 
engines using two alternative analytical methods: the Schmidt analysis and the adi¬ 
abatic analysis. These are two very idealized models yielding far better results than 
those actually attainable in a real world application. Even if the two models are 
very similar, the fact that the Schmidt analysis adopts the additional simplification 
of isothermal expansion and compression spaces induces relevant differences in the 
optimization outcomes, as shown later. 


10.3.2 Objectives of the Optimization 

In this case, the choice of the objectives of the optimization is straightforward. 
In fact, the output parameters we can compute in Schmidt and adiabatic analyses 
are a few general informations over the cycle, like: 

• maximum and minimum temperatures in the cycle, and temperature swing in the 
expansion and the compression spaces, 

• maximum and minimum pressure, and pressure swing in the cycle, 

• pressure phase angle with respect to the expansion volume variation, 

• heat exchanged by the engine, 

• work produced by the engine. 
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Experiments 


Simulations 



1 Objective: 
efficiency 


1 Objective: 
power 


2 Objectives: 
power and 
efficiency 



Direct H Direct H Indirect 

nondimensional I dimensional H dimensional 
parameterization I parameterization I parameterization 



Computation setup: step size and stopping 

criterion 


O pti m i z a ti o n con strain ts 



Fig. 10.6 Summary of the choices made in the setting up of the Stirling engine optimization problem 
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Thus, the choice is quite obvious since the most interesting objectives in engines 
design are the work produced per cycle (or the net power output by multiplying it by 
the revolution speed), and the engine thermodynamic efficiency which is the ratio 
between the net work output and the heat input over one cycle. Schmidt and adiabatic 
analyses actually do not need a revolution speed to be defined. However, we define a 
revolution speed, which is kept constant throughout the whole optimization process, 
and which is just a multiplying factor allowing us to refer to the more commonly 
used power output, in place of the net work output per cycle as output parameter. 

Since the Schmidt analysis always yields the Carnot thermodynamic efficiency, 
in this case we address a single objective optimization aiming at the maximization of 
the net power output. In the case of adiabatic analysis, instead, we address a multi¬ 
objective optimization aiming at the maximization of the net power output and the 
maximization of the thermodynamic efficiency of the engine. 

Also the temperature and the pressure values within the engine are of interest to the 
designer, to avoid excessive thermal and mechanical loads on the engine components. 
For this reason, suitable constraints will be required over these output parameters. 

However, Schmidt and adiabatic analyses just consider the thermodynamic of the 
engine. Other important issues in engines design are not investigated, such as, the 
weight of the engine components, the mechanical stresses, or the size of the heat 
exchangers. 


10.3.3 Input Variables 

The parameterizations required by the Schmidt and the adiabatic analyses are almost 
the same, but Schmidt’s results do not depend on the gas specific heat coefficients, 
while adiabatic outputs do. 

Fourteen parameters are included in the analysis, they are: 

• the swept volumes of the expansion and compression spaces (V sw ,e and V sw ,c)> 

• the dead volumes of the expansion and compression spaces (Vd, e and V^ jC ), 

• the heater, regenerator, and cooler volumes (V^, V r , and \4), 

• the heater and cooler temperatures (Th and 7^), 

• the expansion space to compression space volume phase angle (a), 

• the overall mass of operating fluid (M). This can be substituted by some other 
quantity defining the amount of fluid inside the engine such as, the average cycle 
pressure (p m ), 

• the properties of the operating fluid (R and c p ), 

• the revolution speed ( rpm ). 

The parameterization is more rigid and does not allow many alternative formulations 
as it was in the examples discussed in Chaps. 8 and 9. These parameters can be 
given directly as input or can be given indirectly by defining, for instance, the type 
of Stirling engine, the pistons diameter and stroke, and by computing subsequently 
the swept volumes. Some nondimensional parameters could also be used, however, 
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all the alternatives are absolutely equivalent in terms of simulations and outcome 
of the optimization. The operating fluid which is employed in Stirling engines is 
commonly air, helium, or hydrogen. In general, the larger is the gas constant, the 
better is the engine performance, since a small change in the hot temperature is 
reflected in an elevated pressure driving the pistons motion. Thus, hydrogen is the 
best but it also brings containment problems. Second comes helium which is often 
employed in real engines. Air, despite its relatively low thermodynamic constant, is 
also often used because it is found much more easily in nature, and this makes the 
engine replenishment in case of pressure drop in the working space due to leakages, 
extremely easy. 

We choose to keep the cooler temperature (7^) and the revolution speed ( rpm ) 
constant throughout the optimization process and to consider helium as the operating 
fluid (thus, fixing R and c p ). The remaining parameters are adopted as input variables 
of the optimization. 

We expect that, according to both Schmidt and adiabatic analyses, the optimal 
configurations which will be found when pursuing the maximization of the power 
output, will have approximately zero dead volumes (Vd, e = Vd,c — Vh — V r = 
Vk = 0 cm 3 ) so that the meaningful input variables will actually become five (V sw , e , 
V sw ,o a ’ Th> and Pm )• In fact, both the Schmidt and the adiabatic analyses consider 
isothermal heat exchangers and no constraint is imposed to the heat transfer rate. As 
a consequence, the heat exchanger volumes act as dead volumes to all intents and 
purposes. We also expect the hot temperature and the mean pressure to be as high as 
possible compatibly with the optimization constraints. 

No other parameter except the step size and the stopping criterion for the adiabatic 
analysis needs to be defined for the setup of the simulation process. 


10.3.4 Constraints 

The constraints applied to the Stirling engine optimization problem go beyond the 
typical x m i n <x< x max type. 

Of course we define ranges for the input variables, in particular, we impose: 

• 0cm 3 < V SWie < 400cm 3 , 

• 0cm 3 < V sw ,c < 400cm 3 , 

• 0 cm 3 < Vd, e < 100 cm 3 , 

• 0 cm 3 < Vd, c — 100 cm 3 , 

• 0cm 3 < V/j < 100cm 3 , 

• 0 cm 3 < V r < 100 cm 3 , 

• 0 cm 3 < V k < 100 cm 3 , 

• 350 K < T h < 900 K, 

• T k = 300 K, 

• —7T < a < 7T, 

• 1 bar < p m < 49 bar, 
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• R = 2077 


J 

kg,K’ 


• c p — 5193 kg J K , 

• rpm = 600 rpm. 


Additional constraints are given in order to limit the engine’s size, the stress due to 
the pressure, and the engine minimum power output: 

• V SWf e + Vd,e + Vh + V r + V c + Vd,c + Vsw,c — 500 Cm 3 , 

• Pmax S 50 bar, 

• Pnet > 300 W. 


The last constraint was added to prevent the multi-objective optimization to move 
towards zero power output configurations, which are likely to happen pursuing the 
objective of maximum thermodynamic efficiency. 

The ranges of the input variables are restricted as the optimization process goes on. 


10.3.5 The Chosen Optimization Process 

A similar optimization process is applied twice, using the Schmidt analysis first, 
and then using the adiabatic analysis. The simulation process is a cheap analytical 
computation which requires a fraction of a second to be completed on a personal 
computer. For this reason, this optimization exercise is also used for comparing 
different optimization methods. 

We start from considering the Schmidt analysis for the ten input variables problem. 
At first a Sobol DOE with 2048 feasible designs is performed. The range of the 
input variables is then restricted around the optimum solution found, and a 1P1-ES 
stochastic optimization with 1024 designs is applied. The range of the input variables 
is restricted once again around the optimum configuration found, and a Nelder and 
Mead simplex deterministic algorithm is applied in the end. 

As for the adiabatic analysis for the ten input variables problem, A Sobol DOE 
with 2048 feasible designs is performed first. Then the range of the input variables 
is restricted, and a MOGA with 4096 designs (32 individuals x 128 generations) is 
applied. The MOGA is followed by two 1P1-ES with 1024 designs each: the first 
aiming at the maximization of the power output, the second aiming at the maxi¬ 
mization of the engine’s thermodynamic efficiency. The two 1P1-ES are followed 
by two Nelder and Mead simplex optimizations having the same objectives as the 
evolutionary optimizations. 

Thus, the procedures followed for the Schmidt and the adiabatic analyses are 
much the same. The differences are that: 

• no MOGA is performed using Schmidt analysis, since for the Schmidt case the 
thermodynamic efficiency objective loses its significance, 

• the evolutionary and the simplex steps are performed twice, once for each opti¬ 
mization objective, in the adiabatic analysis. 




10.3 Methodological Aspects 
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As expected, the results of the optimization tend to lead to configurations with zero 
dead volume, maximum heater temperature, and maximum total volume, where the 
total volume is 


Vtot — Vsw,e + Vd,e + Vh + V r + Vk + Vd,c T Vsw,c- (10.35) 

For this reason, in the second part of the optimization process, the heater tempera¬ 
ture is fixed to 900 K, the dead volumes to zero, and the compression swept space 
to V sw ,c — 500 cm 3 — V S w,e- In this way, we define an optimization problem whose 
three variables are: V SWte , a, Pmean- Over the new optimization problem the same 
optimization process adopted for the ten input variables case is applied, thus involv¬ 
ing: a Sobol DOE, a MOGA optimization, a 1P1-ES optimization, and a Nelder 
and Mead simplex optimization. Each optimization algorithm was initialized from 
the best configurations found in the previous step of the process. Actually, the zero 
dead volume condition is approached only when the net power output objective is 
addressed. When the thermodynamic efficiency is addressed, the optimal solutions 
present large dead volumes and very poor performance in terms of power output. 
The reason for this, is that the thermodynamic efficiency reduction in the adiabatic 
analysis is due to the non-isothermal behaviour of the expansion and the compres¬ 
sion spaces. Temperature variations in the expansion and the compression spaces 
are mainly due to the pressure variation in the working space caused by the pistons 
motion. For this reason, the thermodynamic efficiency is high when 

• the swept volumes are low, since a low swept volume means also a low pressure 
variation in the working space over the cycle, 

• the dead volumes are large, since the dead volumes act as a buffer volume con¬ 
taining the pressure and temperature variations in the working space. 

These two conditions heavily and negatively affect the engine performance in terms 
of power output. Since the simulations require a very short computing time to be 
completed, in the last part of the optimization process, several optimization tech¬ 
niques were compared starting from scratch, using the same initial design point or 
population, and the same design space. Single objectives techniques were compared 
over the Schmidt analysis, and multi-objective techniques were compared over the 
adiabatic analysis. The comparison also included a few DOE+RSM techniques cou¬ 
pled to metamodel-assisted optimization processes. 

A summary of the elements involved in the first two parts of the optimization 
process is given in Fig. 10.7. 


10.4 Results 

Let us consider the ten variables Schmidt optimization problem. The results of the 
Sobol DOE, the 1P1 -ES, and the simplex optimizations are summarized in Table 10.1. 
The whole optimization process is carried out by using the optimization dedicated 
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Fig. 10.7 Elements involved in the Stirling engine optimization problem 


Table 10.1 Ten input variables optimization: Schmidt analysis results 


Input Variable Sobol DOE (2048 designs) 1P1-ES (1024 designs) Simplex (1083 designs) 
or Output Range Best Range Best Range Best 

Parameter Low High P net Low High P net Low High P net 


a [deg] 

-180 

180 

123.53 

10 

170 

110.05 

80 

150 

111.91 

T h [K] 

350 

900 

839.51 

400 

900 

888.39 

600 

900 

899.96 

Vsw,e [cm 3 ] 

0 

400 

288.62 

10 

350 

319.86 

50 

350 

329.92 

V sw ,c [cm 3 ] 

0 

400 

114.40 

10 

350 

175.04 

50 

350 

169.41 

Vd, e [cm 3 ] 

0 

100 

22.60 

0 

100 

1.73 

0 

50 

0.01 

Vd,c [cm 3 ] 

0 

100 

8.58 

0 

100 

0.00 

0 

50 

0.23 

V h [cm 3 ] 

0 

100 

24.60 

0 

100 

2.66 

0 

50 

0.11 

V r [cm 3 ] 

0 

100 

21.35 

0 

100 

0.55 

0 

50 

0.30 

V k [cm 3 ] 

0 

100 

5.60 

0 

100 

0.00 

0 

50 

0.00 

Pm [bar] 

1 

49 

29.26 

4 

46 

24.94 

15 

40 

25.58 

M[ g] 

— 

— 

0.671 

— 

— 

0.460 

— 

— 

0.466 

Pmax [bar] 

— 

— 

41.69 

— 

— 

49.70 

— 

— 

50.00 

Vtot [cm 3 ] 

— 

— 

485.75 

— 

— 

499.85 

— 

— 

499.99 

Pnet [kW] 

— 

— 

2.756 

— 

— 

5.285 

— 

— 

5.474 

7?[%] 

— 

— 

64.26 

— 

— 

66.23 

— 

— 

66.67 


software modeFRONTlER. At each step, the design space size is shrinked around 
the best configuration found in the previous step. As expected, the optimization is 
clearly moving towards an optimum configuration with zero dead volumes, heater 
temperature of 900 K, total volume of 500 cm 3 , and maximum pressure in the cycle 
of 50 bar. Despite 2048 configurations were evaluated in the Sobol DOE, the best 
result found by the process is still far from the optimality condition. In fact, the 
number of input variables is rather large, and a deep investigation of the design space 
is not attained even with such a number of simulations. 

A stochastic optimization is more precise in finding optimum solutions than 
pseudo-random searches. In fact, passing from the Sobol sampling to the 1P1-ES 
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the performance of the better configuration in terms of maximum net power output is 
almost double. Deterministic optimization is even more precise than stochastic opti¬ 
mization, and the performance of the better configuration is further improved after 
the simplex optimization. Thus, the procedure starts from a quasi-random exploration 
of the design space and moves on, step by step, towards an accurate refinement of 
the solution. 

The choice of any optimization process is always a trade-off between how much 
importance is given to the design space exploration and to the solution refinement, 
that is, between robustness and velocity. By robustness the capability of avoiding 
local optima and explore the whole design space is meant. 

The same procedure is followed for the ten variables adiabatic optimization prob¬ 
lem, whose results are summarized in Table 10.2. Since in the adiabatic analysis 
we address two objectives, we also apply multi-objective optimization algorithms. 
The results in terms of maximization of the net power output (right hand of the 
Pareto frontiers in Fig. 10.8) go in the same direction of those found for the Schmidt 
analysis (zero dead volume, maximum heater temperature, maximum total volume, 
maximum average pressure compatibly with the maximum pressure constraint). As 
for the maximum thermodynamic efficiency objective, as already noted, it must be 
considered that the source of inefficiency in the adiabatic analysis is due to the non- 
isothermal behaviour of the expansion and the compression spaces. Thus, the least 
is the temperature variation in those spaces, the better is the efficiency. However, 
the least is the temperature variation, the worse is the net power output due to the 
small working space compression needed for causing small temperature variations. 
An additional constraint on the minimum net power required from the engine is 
given in order to avoid degenerate solutions. Despite this constraint, it is clear from 
Table 10.2 that the optimum solutions in terms of thermodynamic efficiency have 
small compression and expansion space swept volume, elevated dead volume, high 
mean pressure, low pressure and temperatures swing in the working space over the 
cycle. In other words, if it was not for the constraint given on the net power output, 
the optimum configuration would have moved towards an engine which stands still 
and, obviously, gives no power output. 

In the second part of the optimization procedure 

• the dead volumes are fixed to zero, 

• the total volume is fixed to 500 cm 3 , 

• the heater temperature is fixed to 900 K, 

that is, the number of input variables is reduced to three. Actually, the dead vol¬ 
umes are fixed to 1 mm 3 each, since zero dead volumes causes the adiabatic analysis 
to diverge. The design space is further reduced at each step in a neighbourhood of 
the optimum solutions previously found. Table 10.3 shows the results of the three 
input variables Schmidt optimization; Table 10.4 shows the results of the three input 
variables adiabatic optimization. 

The adiabatic optimization has a larger design space since it must follow the ten¬ 
dencies of the two objectives of the optimization. Now that there is no dead volume 
to play with, when pursuing the maximization of the thermodynamic efficiency, the 
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10.4 Results 


Fig. 10.8 Evolution of 
the Pareto frontier for 
the Stirling engine adi¬ 
abatic analysis MOGA 
optimization 



results of the optimization find a different strategy for limiting the temperature swing 
in the expansion and the compression spaces. This strategy tends to promote high 
values of a , as demonstrated by the results in Table 10.4. In fact, a high value of a 
means that the volume variations in the expansion and the compression spaces are 
almost in counterphase so that when the expansion space is large, the compression 
space is small and viceversa. Overall, the size of the working space (V e (0) + V c ( 9 )) 
is not undergoing large variations over the cycle. As a result the ratios and 

are reduced. Simplex algorithm brought no improvement at all for the adia- 

‘ m i n 

batic analysis optimization, and a very small improvement for the Schmidt analysis 
optimization. This could mean that the 1P1-ES algorithm had already reached the 
optimum solution, at least locally. 

The results of the comparison between different single objective algorithms over 
the three variables Schmidt optimization are shown if Fig. 10.9 and Table 10.5. All 
the algorithms were started from the same initial point. 

BFGS and Fevenberg-Marquardt algorithms fail to converge to the optimum 
solution. The reason for this failure is the same already discussed in the context of 
Example 4.1. In fact, the initial point of the optimization is near to the border of the 
feasible region, since the value of its maximum pressure over the cycle is almost 
50 bar. However, the gradient pushes the algorithm to increase the mean pressure in 
the cycle because the mean pressure is proportional to the net power output. In this 
way, the maximum pressure constraint is violated, the objective function is penal¬ 
ized, and the gradient estimation is incorrect. As a result, the algorithms get stuck 
almost immediately. This shows that BFGS and Fevenberg-Marquardt algorithms, 
even being very effective, only work properly in unconstrained optimization; their 
application to constrained optimization problems is likely to fail as soon as a con¬ 
straint is violated during the optimization process. The remaining algorithms have 
comparable efficiency. 1P1-ES shows a slower convergence rate, while the DES is 
not only faster than 1P1-ES, as expected, but, surprisingly, is also almost as fast 
as the simplex and the NFPQFP algorithms. NFPQFP encounters problems in the 
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Iterations 

Fig. 10.9 Convergence speed of different methods over the three variables Stirling engine opti¬ 
mization problem through Schmidt analysis 


first iterations and still yields no improvement after 35 iterations, but later it quickly 
makes up for the lost time and at iteration 70 is leading over all other algorithms. 

Trying to generalize the above observations, let us suppose the simulations we are 
running are computationally intensive and each iteration requires 6 h of CPU time 
to complete. The deterministic optimization algorithms, unless they were failing, 
would have required between two and three weeks time (56-84 simulations) to reach 
a reasonably good approximation of the optimum configuration, and at least one 
month (120 simulations) to meet the stopping criterion and terminate. If we want 
to speed up the optimization process we could barter accuracy for speed using a 
DOE+RSM approach. We have tried two different DOE+RSM approaches: 

• a 2-levels full factorial DOE (8 designs) plus a 3 levels central composite faced 
DOE (7 designs) coupled to a Gaussian process response surface, 

• a uniformly distributed latin hypercube using 32 designs coupled to a Kriging 
response surface. 

According to our hypothesis the first would have required less than 4 days of CPU 
time (let us say, a long week end, from friday late afternoon to tuesday morning), and 
the second 8 days. An optimization process could then be applied to the metamodel 
running in a fraction of a second. Fortunately the Schmidt analysis requires less 
than 0.1 sec to complete and we do not have to worry about CPU time, actually. 
However, the two DOE+RSM processes discussed above gave amazingly accurate 
results considering the small number of simuations they required (see Fig. 10.9 and 
Table 10.5). 

It is true that using metamodels means to accept some potential degree of inaccu¬ 
racy in predictions, however, it can also save a lot of time sometimes. The results of 
the comparison between different multi-objective algorithms over the three variables 
adiabatic optimization are shown in Table 10.6 and Fig. 10.10. 
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Fig. 10.10 Pareto frontier evolution for different methods over the three variables Stirling engine 
optimization problem through adiabatic analysis 


MOGA, MOSA, and NSGA were started from the same population of 12 indi¬ 
vidual obtained with a Sobol DOE and ran for 84 generations so that the number 
of simulations was reaching 1, 000 at the end of the optimization. Note that a large 
part of the design space defined by the ranges in Table 10.6 is considered “unfea¬ 
sible” since it causes the results of the adiabatic analysis to break the maximum 
pressure constraint. For instance, out of the 12 individuals of the initial population 
10 were unfeasible. The MMES was started from a population of 4 individuals (2 
feasible and 2 unfeasible) taken from the MOGA initial population, and ran for 50 
generations using an adaptive (4, 20)-ES scheme with maximum individual life span 
of 5 iterations. MOGT was started from a feasible individual of the MOGA initial 
population. The remaining parameters for the setup of the optimization were as in 
Example 5.1, except for the MOSA which was started from a temperature of 0.2 and 
had ^ fraction of hot iterations. 

Due to the strong limitations caused by the maximum pressure constraint, MOGT 
failed to converge and stopped after 62 simulations, 59 of which were unfeasible. 
The incidence of unfeasible samples is limited to 18 % in MOGA, 22 % in MOSA, 
25 % in NSGA, and 9.5 % in MMES. Apart from MOGT, the other algorithms show 
a rather good convergence towards the Pareto frontier (see Fig. 10.10), MMES being 
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(a) nlpqlp optimum configuration after Tab. 10.5: 
phase plane volumes and pressure plots 
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(b) NSGA max Pnet optimum configuration after (c) NSGA max r] optimum configuration after Tab. 10.6: 
Tab. 10.6: phase plane volumes, pressure, and tem- phase plane volumes, pressure, and temperatures plots 
peratures plots 




(d) NSGA max Pnet optimum configuration after (e) NSGA max 77 optimum configuration after Tab. 10.6: 
Tab. 10.6: work and heat quantities exchanged over work and heat quantities exchanged over one cycle 
one cycle 


Fig. 10.11 Thermodynamic cycle details for some optimum configurations after Tables 10.5 and 
10.6. The x signs in the phase plane plots are placed every A 0 = 90° crank angle, the arrow 
individuates the 0 = 0 ° locations and the direction in which the path is travelled 


slightly less efficient than the other methods towards the end of the process (after 
iteration 500), and MOGA being less efficient at the beginning of the process (before 
iteration 200). The final Pareto population is made of 43 individuals for MOGA, 27 
individuals for MOSA, 66 individuals for NSGA, and 49 individuals for MMES. 
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Figure 10.11 shows some information about the thermodynamic cycle of the 
optimum configuration reported in Table 10.5 found by the NLPQLP optimization 
algorithm through Schmidt analysis, and by the optimum configurations reported in 
Table 10.6 found by the NSGA optimization algorithm through adiabatic analysis. 

The irregular shape of the temperature phase plot in Fig. 10.1 lb, c is due to the fact 
that the expansion and the compression space dead volumes are zero. The process is 
as follows: we consider the expansion space, as initially containing a certain amount 
of fluid at a certain temperature; as the space is reduced to zero the fluid is ejected 
from the space completely; fresh fluid enters the space when the piston recedes. The 
incoming fluid, however, has temperature 7^ due to the conditional temperature in 
Eq. 10.18, no matter what the temperature of the operating fluid in the space was 
before. This originates the discontinuities in the fluid temperature observed in the 
expansion space, and in the compression space. The effect would have been avoided 
if the dead volumes were present, since dead volumes act as buffers, hence smoothing 
out the sudden temperature change in the spaces due to the incoming operating fluid 
from the heat exchangers. 

Note that, if the adiabatic simulation has reached convergence the energy balance 
equations are fulfilled: at the end of the cycle we have j> Q r = 0, § Qh = § W e , and 
§ Qk = § W c (see Fig. 10.1 Id, e). 


10.5 Conclusions 

The optimization of Stirling engines has been addressed by using the Schmidt and the 
adiabatic analyses. The optimization process was performed for each type of analysis. 
The process was quite standard and involved a Sobol DOE, a MOGA, a 1P1-ES, 
and a simplex optimization. After the optimization process was completed, it was 
noted that some input variables were moving towards one extremity of their range. 
When such a behaviour is found, it is clear that the input variables would move even 
furhter if they were not constrained by their ranges. Under these circumstances, two 
possible choices are suggested: 

• if possible, move the ranges to comply with the tendencies of the input variables, 
this could lead to better performing solutions, 

• if not possible, change the input variables to constants and proceed with the opti¬ 
mization process. 

In our case, since negative volumes have no physical meaning, and a higher heater 
temperature would have damaged the engine, we can not move the variables ranges 
and we choose the second possibility. In this way, we help the optimizer to man- 
tain the optimum values for some of the variables. In fact, due to the randomness 
which is present in stochastic optimization the optimizer found difficulties in keeping 
the values of these variables anchored to the extremity, thus wasting time running 
suboptimal simulations. Moreover we can now proceed with an easier optimization 
task which will run faster since it involves a lower number of input variables. After 
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the number of input variables was decreased, the same optimization process was 
repeated for the “reduced” problem. This strategy was successful in reaching an 
optimum solution and to further improve it. This latter task was achieved by both 
reducing the number of variables, and by progressively moving from an initial explo¬ 
ration of the design space by means of a Sobol DOE or a MOGA optimization, to a 
refinement of the solution through a Nelder and Mead simplex optimization. 

The results given by the two analyses in terms of optimum engines are quite 
different from each other, not only for the fact that the Schmidt optimization is 
single objective and the adiabatic optimization is two-objectives. For instance, let 
us consider the optimum configurations in terms of maximization of the net power 
output from Tables 10.3 and 10.4 

• for the Schmidt analysis we have: 

a = 113.5°, V SWy e = 342.8cm 3 , V sw , c = 157.2cm 3 , p m = 26.4bar, P net = 
5.513 kW, 

• for the adiabatic analysis we have: 

a = 146.2°, V sw<e = 329.5cm 3 , V sw<c = 170.5cm 3 , p m = 30.2bar, P net = 
4.274 kW. 

The main difference in the two configurations is in the larger a value which is attained 
in the adiabatic analysis. This causes the compression rate to be reduced over 
the cycle, which allows higher mean pressures to be applied to the engine without 
breaking the maximum pressure constraint. It also causes the ellipse in the volumes 
phase plot in Fig. 10.11 to be more elongated and tilted towards the left side of the 
plot. 

From the point of view of optimization it is interesting the comparison between 
different methods, as discussed in the previous section and summarized in Tables 10.5 
and 10.6, and in Figs. 10.9, 10.10, and 10.11. 

It could be argued that too many constraints were applied to the optimization 
problem. For instance, the P net > 0.3 kW was somewhat limiting the maximum 
efficiency solutions which were found. However, as it was already reminded before, 
this constraint was due to the fact that it was clear that the engine would have moved 
towards a degenerate solution. The a < 160 ° constraint applied in the end was lim¬ 
iting even more the action of the optimization algorithms in terms of maximization 
of the thermodynamic efficiency. In fact, all the maximum efficiency optimum solu¬ 
tions in Tables 10.4 and 10.6 have the value of a at 160 0 or very close to 160 °. This 
indicates that the optimizer would have gone further if it could, tilting and elongating 
even more the volume phase plot in Fig. 10.11, increasing p m towards 50 bar, and 
reducing the net power output up to when the P net >0.3 kW constraint was met. 

It would be interesting to investigate the Stirling engine optimization problem 
using more advanced and realistic simulation models. The results found in this exer¬ 
cise for a are quite similar to the values which are used in real engines, while the 
swept volumes ratio *4^, in real engines, is generally not too far from 1. In our 

v sw,c 

analysis, instead, > 2 was found. 

v sw,c 





